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Abstract 

A' pseudospectral explicit scheme for solving linear, periodic, parabolic 
problems is described. It has infinite accuracy both in time and in space. 
The high accuracy is achieved while the time resolution parameter M 
(M — 0( ^ ) for time marching algorithm) and the space resolution parameter N 
(N = 0( s )) have to satisfy M = 0(N* +e ) e > 0, compared to the common 
stability condition M = 0(N 2 ) which has to be satisfied in any explicit 
finite order time algorthm. 
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Introduction. 


In recent years, it has been shown that spectral methods can provide 
a very useful tool for the solution of time dependent partial differential 
equations [3], A standard scheme uses spectral methods to approximate 
the space derivatives and a finite difference approach to march the 
solution in time. This tactic results in an unbalanced scheme; it has 
infinite accuracy in space and finite accuracy in time. It is obvious 
that the overall accuracy is influenced strongly by the relatively poor 
approximation of the time derivative. Moreover, using finite order 
explicit scheme results in a very stringent stability condition. The time 
step. At, has to satisfy 

At = 0(i=-) (1.1) 

N 

where N is the number of grid points in space. This severe condition 
is commonly overcome by resorting to implicit schemes. Varga [6], Cody, 
Meinardus and Varga [2] approached these problems by using Chebychev 
rational approximations of the evolution operator. Thus, they overcome 
two drawbacks - low accuracy and stringent stability condition. In 
fact, the implicit scheme presented in [2], [6] is unconditionally stable, 
and the error in time decays exponentially. 

Implicit algorithms involve inverting matrices. When the space 
approximation is based on finite differences or finite elements (as in 
[ 2 ] * [ 6 ])> the related matrices are banded ones (e.g. tridiagonal) which 
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wakes them relatively easy to invert. On the other hand, using spectral 
methods for the space discretization results in full matrices. Inverting 
these matrices is a time consuming procedure. 

In this paper we describe an explicit scheme for the solution of 
parabolic problems when the space discretization is done by spectral 
methods. This scheme is highly efficient (its efficiency is equivalent 
to having a stability condition At = O(^j-) ) and the error in time decays 
exponentially. In Section 2 we present a model problem and its fully 
discrete solution. The new approach for approximating the evolution 
operator is described in Section 3. In Section 4 we carry out an error 
and stability analysis. Numerical experiments confirming the theoretical 
results are presented in Section 5. 
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2. The Model Problem. 


Let us consider the heat equation 


U - GU = 0 


o < x < n 


U(x, 0) = IT 00 


( 2 . 1 ) 


u(o, t) = u(n, t) = o 


where G is the spatial operator 


G = a- 


( 2 . 2 ) 


Discretizing (2.1) in space using pseudospectral Fourier method results 
in a semidiscrete representation 


( V - G N U N = 0 


U N (x, 0) = U“(x) 


(2.3) 


u N (0, t) = u N (n, t) = o 


U N =P N U > G N =P N GP N ! U N =P N U ° 


(2.4) 


and where for any function f(x), P^f(x) is its sine interpolant at the 
collocation points 


Xj = jn/N 


j = 0, 1,..., N-l 


(2.5) 
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or more precisely. 


N-l 

P M f(x) = l a, sin(kx) (2.6) 

N k=0 K 

where 

, N-l 

\ = M I f(x )sin(kx.) . (2.7) 

K N j=0 J 3 


G^ is an operator defined on N dimensional subspace ; thus it can 
be represented as a N * N matrix. The formal solution of (2.3) is 


U N (x,t) = exp(tG N )ujj(x) 


( 2 . 8 ) 


where exp(tG^) is the exact evolution operator. A fully discrete solution 
of (2.1) is achieved by approximating this evolution operator. In 
[5] , it has been shown that any explicit time scheme can be represented 
as 

^ = H M ttG lX ' 2 ' 9 > 


where H M (z) is a polynomial of degree M which converges to e in 
the domain which includes all the eigenvalues of the operator tG^ . 
is the fully discrete solution and H^(tG^) is the numerical 


evolution operator. 
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3. The Orthogonal Polynomials Scheme. 


Let E be the error that results from approximating the evolution 


operator. Then 


E = [exp(tG N ) - H M (tG N )]U^ . 


(3.1) 


The eigenvectors of the matrix tG^ are W^,...,W N where 

(W,) = sinfkx.) . Due to the orthogonality of this set of eigenvectors, 

j ’ 

N is a normal matrix and there is an orthogonal matrix S^ T such that 


e = s n s"V. 

N N N N 


(3.2) 


where D is the diagonal matrix 


<V kk ■ e ‘ -w 


(3.3) 


and A^t are the eigenvalues of tG^. Since is orthogonal matrix, 


we have S* 


S’ |=1. Therefore, 


|E| l L S ||S N || M D N || | |S N *| | = 
2 L 2 L 2 L 2 


I E | |, * max |e - H (z) 

L 2 z€I 


(3.4) 


where I is the domain which includes all the eigenvalues of tG^. 


For our case 


I = [-aN 2 t,0] . 


(3.5) 
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A standard finite order scheme can be characterized by a polynomial 

z 

H^(z) based on a Taylor expansion of e . Thus, it has high accuracy 
only for a small z . The error increases rapidly when z is increased. 
This property explains the poor accuracy and stringent stability condition 
mentioned in the introduction. 

Let us take, for example, the Modified Euler scheme. The numerical 
evolution operator is 

H M (t V " 0 * AtG N * 1 (4t S )2 ) f3 ' 6) 

where 

At = t/n . (3.7) 


Thus 



Figure 1 


(3.10) 
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From Figure 1, we find that H,,(z) converges to e Z when 
-2n £ z <: 0 . 


(For accuracy, a more stringent condition is necessary). Using (3.5), 
(3.7), and (3.10) results in the following stability condition 

At * \ • (3.11) 

N 

Expression (3.4) suggests that a uniform approximation of e Z is 
preferable. Such an approximation is achieved when one uses Chebychev 
polynomials expansion of the exponential function (see discussion in [5] 
for hyperbolic problems) . Let 

w = i(z + R) -UwM (3.12) 

where 

R = j aN 2 t . (3.13) 

It then follows 

e Z = e" R e Rw = £ b T (w) (3.14) 

k=0 K 

where T k (w) is the Chebychev polynomial of order k and [1] 

1 - 1/2 

b k = e ~\ / e Rw T k (w) (1 - w 2 ) dw = e' R c k I k (R) 


(3.15) 



and also 


c k ~ 


k = 0 


k 5 1 


(3.16) 


I, (R) is the modified Bessel function of order k. 
k 

polynomial approximation of e Z is 


Thus, the M degree 


M 

H M (z) = ^ W w(z)) • 

M k=0 k 


(3.17) 


Since (3.12) we substitute the operator defined as 


f n = »!"»* RI ’ 


(3.18) 


for w. H^(F^) is the numerical evolution operator. Thus, the fully 


discrete numerical solution of (2.1) is 

M 


V N ' W U S * l WV U S 

k=0 


(3.19) 


?k ( Ffj) is com P uted by using the recurrence relation 


Hence 


T k (x) = 2xT k _ 1 (x) - T k _ 2 (x) 


T Q (x) = 1 ; Tj(x) = x . 


k £ 2 


T k<V U S - 2F N T k-l"V U N - T k-2t F N )U S * 4 2 


(3.20) 


V F N> U S ■ U N 


■V f n )u n ' f n u n° 


(3.21) 
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The algorithm defined by (3.19), (3.21) can be regarded as a three 
level scheme since it uses the recurrence relation. Therefore, it has 
the disadvantage of requiring extra memory. There are two possible 
ways to overcome this drawback. The first one is to convert (3.19) to 
a power series in F N and using Horner scheme to compute vj^. The 
disadvantage of this approach is its sensitivity to round-off errors. 
The second one is based on calculating the roots of H^(w) . Let us 
assume that the roots are 


© ^ 9 * * * 9 * (3.22) 

Since the b^ are real, every complex root appears with its conjugate. 
Rearranging (3. 22| is such a way that the first 2p roots are p conju- 
gate pairs, we get 
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Hence we get 


w = 


p 

n 

i=l 


[I - o.F m 
1 l N 


Vn> 


M-p 

n 

i=2p+l 


[I - V t F N l 


uj; (3.26) 


Each one of the algorithms described above can be used as a one 
step method by calculating the solution at the final time t directly 
from the initial data. It can also be used as a marching scheme if one 
is interested in intermediate results. The size of the time step At 
depends only on the information one wants to get out of the numerical 
procedure. At enters instead of t in the expressions above, and the 
parameter R is determined accordingly. In any case, the refinement of 
the algorithm is done by increasing the degree of the polynomial and not 
by decreasing the size of the time step. 
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4. Accuracy and Stability. 


Using (3.4), (3.15), and (3.17) we get 

oo 

•Ri 


|E| | S 2e 
L 2 


l I (R)T (w) 
k=M+l K K 


-H w J 1. 


(4.1) 


Rw 

Since e is an entire function it satisfies the following theorem 
([4], pp. 94-96). 


Theorem . (S . N. Bernstein): Let f(w) be an entire transcendental 
function which is real for real w . Then there exists a sequence of 
integers n^, n^,... with n^ 00 such that the relation 


n v ' 

lim = 1 

l« n + j| 

V 


(4.2) 


holds, where are t ^ e coefficients in the expansion 


f(w) = + l a k T k (w) 

k=l 


(4.3) 


and 


E n (f) = 


f(w) - - l CX, T (w) 

L k=l k k 


(4.4) 


There is a sequence of integers n^ , V = 1, 2,... of the above 
type provided 

1. °n +1 ^ 0 v = 1* 2,... (4.5) 

V 


and 
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2. I |a | = 0 (| a |) as y “ 

k=n +2 k V 

y 


(4 


In our case we can take n^ = y , y = 1, 2,.. 


and it follows that 


I e Mi * 2e"\. .(R) (1 + 0(1)) 


M+l 


(4. 


The asymptotic expansion of I^(R) is [1] 


R 


1,00 --1= U . , (EOHy-9) 


sn . r 


2! (8R) 


(y-1) Cy-9) (y-25) 
3!(8R) 3 


where 


y = 4k . 


(4, 


(4 . 8) can be written as 


2e_R v r > 


i - 


iL_ + JL (U_-) 

8R 2! l 8R J 


. + 0 ( 


*] 


(4. 


or 


2e" R I k (R) 'v/^l exp(-y/8R) + O^R ^ 


(4. 


From (4.7)-, (4.9) and (4.11) we conclude that an e time accuracy. 


E|L Se , 
L 2 


(4. 


is achieved when 


M = 0(R 1/2 ) . 


. 6 ) 


7) 


(4.8) 


9) 


10 ) 


11 ) 
c.e . 
12 ) 


(4.13) 
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It is clear that satisfying (4. 13) guarantees stability. In fact, using 
(3. 1) , (4. 12) we get 

||exp(tG N ) -H^ttyll s e ; (4.14) 

hence 

* | | exp(tG N ) | | + e. (4.15) 

Since exp(tG N ) is a stable operator [3], H M (tG N ) is stable as well. 

2 

R is equal to aN t/2 ; thus from (4.13) we can conclude the main 
result of this ^analysis: In order to achieve e time accurate, stable 
solution of (2.3), M has to satisfy 

M = 0(N) . (4.16) 

A similar analysis for any finite order scheme based on Taylor expansion 

of e Z will imply that M [M = O(^) , see (3.7) - (3.8)] has to be 
2 

proportional to N ; thus the advantage of the orthogonal polynomials 
approach is obvious. 


Algorithm Refinement . 

From (3.13), (4.7), (4,9), and (4.11) we get 

2 1 ( 2 \ 

E W N fe 5 exp(-(M/N) Z /at) • 


C4.17) 


Expression (4.17) suggests that refinement of the algorithm while 



14 


M = N° (a > 1) (4.18) 

will yield an exponential decay of the error. The accuracy thus 
achieved is the desired spectral accuracy. 


5. Numerical Results . 

Table 1 presents the stability properties of the 0. P. S. 
(Orthogonal Polynomial Scheme) compared to modified Euler scheme 
which is second order in time. We have used the model problem (2.1) 
with a = 1 , and initial data 

U^(x) = sin(3x) . (5.1) 

The solution is computed at t = 1. 


Table 1 



Modified Euler 

0. P. S. 

N 

M 

M 

16 

48 

24 

32 

192 

48 

64 

768 

96 


N - Number of grid points. 

M - The degree of the evolution operator. 
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M indicates the minimal number of applications of the operator tG^ 
one has to use in order to achieve stable (meaningful) results. 

The second table clarifies the spectral convergence of the O.P.S. scheme. 
(We included in this table the results for the Modified Euler scheme as 
well, for the sake of comparison.) The problem solved is 

U + - U =0 
t xx 

0 $ x $ 2 tt . (5.2) 

l/*(x) = x(x - 2ir) 

Note that the periodic continuation of U°(x) belongs to C°; thus 
the Fourier coefficients of U°(x) are decaying slowly. The solution 
is computed at t = 1. 


Table 2 


1 

N 

— 

Modified Euler 

0. P. S. 


M 

1*2 Error 

Ratio 

M 

L 0 Error 

Ratio 

16 

62 

.3791-04 


26 

.1026-04 





17,4 



92 

32 

250 

.2126-05 


61 

.1107-06 





16.2 



134 

64 

1000 

.1339-06 


140 

.8263-09 
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The refinement of the Modified Euler scheme is done while M satisfies 
M = 0.97 x (N/2) 2 • 

For the 0. P. S. algorithm, M satisfies 
M = 2.5 x (N/2) 1 • 2 . 

The increasing ratio between the L^-errors of two successive refinements 
verifies the spectral convergence of the 0. P. S. algorithm. 

In Table 3 we compare the 0. P. S. to the modified Euler scheme from 
the point of view of the amount of work needed to achieve a certain degree of 
accuracy. The problem solved is U - 1)^ = 0 with U (x) = sin(3x) . 

The L^-Error is computed at the time level t = 1, and the space resolu- 
tion is N = 32. 


Table 3 


L 2 

Error 

M (Modified Euler) 

M (0. P. S.) 

1.3 

10- 2 

200 

50 

1.3 

O 

I 

2000 

60 

1.3 

IQ' 6 

20000 

70 
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Conclusion. 


The algorithm presented in this paper achieves the goal of spectral 
accuracy in time and space for the simple model problem C2.1). We believe 
that this approach can be useful for more complicated problems. In fact 
the scheme described in Section 3 is applicable whenever one can represent 
the solution as exp(tG^)U^ and the eigenvalues of tG^ are grouped 
close to the real axis. 
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